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We calculate the spin relaxation time of mobile electrons due to spin precession between random 
impurity scattering (D'yakonov-Perel' mechanism) in electrically gated bilayer graphene analytically 
and numerically. Due to the trigonal warping of the bandstructure, the spin relaxation time exhibits 
an interesting non-monotonic behavior as a function of both the Fermi energy and the interlayer 
bias potential. Our results are in good agreement with recent four-probe measurements of the 
spin relaxation time in bilayer graphene and indicate the possibility of an electrically-switched spin 
device. 



Many fascinating properties of electrons in graphene 
have been brought to light since its discovery, such as 
their high electron mobility and the emergence of anoma- 
lous integer quantum Hall plateaus [1-3]. One of the 
less studied but important questions is the capability of 
graphene to store and transport electron spin. Compared 
with semiconductors such as Si or the III-V compounds, 
graphene bears superior traits for long spin coherence: 
its low density of nuclear spins reduces hyperfine interac- 
tions that are limiting spin coherence in GaAs, while its 
low atomic weight implies intrinsically weak spin-orbit 
interactions (SOI) thus allowing for slow spin relaxation 
[4]- 

Graphene spin valve devices have been demonstrated 
soon after the discovery of graphene [5], followed by 
four-probe spin transport experiments using ferromag- 
netic cobalt [6] and permalloy [7] electrodes. From Hanlc 
precession measurements, spin relaxation times on the 
order of 150 ps were found, and by simultaneously mod- 
ifying the mobility and spin relaxation time by tuning 
the Fermi energy with an external gate, a behavior of 
the spin relaxation time consistent with an Elliot- Yafet 
type mechanism was identified [6]. Subsequent experi- 
ments have confirmed this for single-layer graphene [8], 
but D'yakonov-Perel' type behavior in combination with 
spin relaxation times up to a few nanoseconds at 4 Kelvin 
were found in bilayer graphene [8, 9]. 

Motivated by these observations, we calculate the spin 
relaxation rate for bilayer graphene according to the 
D'yakonov'-Perel mechanism. Our starting point is the 
band Hamiltonian of AB-stacked bilayer graphene (BLG) 
for momenta Kk. = h(k x ,k y ) near the Dirac points K 
(t = 1) and K' (r = -1) [10], 
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in the basis A\, B2, Bi, A2 where A\ refers to the A- 
sublattice in the upper layer, B2 to the B-sublattice in 
the lower layer, etc., and where p — h(k x +iTk y ) — hke lT ^ 
with = a,rct&n(k y / k x ) . Here, the intralayer hopping 
parameter 70 = 2.8 eV determines the Fermi velocity 



vp = ^ajo/h = 8.0 • 10 m/s, whereas the interlayer hop- 
ping parameter 71 = 0.39 eV gives rise to a strong cou- 
pling of the two stacked lattice sites B\ and A2. Skew 
interlayer hopping with strength 73 = 0.315 eV intro- 
duces an additional velocity 1*3 = \a^2,jh, = 5.9 • 10 4 m/s 
and causes a significant trigonal warping of the energy 
dispersion. A tunable energy offset U between the two 
layers can be achieved by applying a bias voltage and 
leads to the opening of a band gap, which has been ob- 
served to reach up to 250 meV [11]. For what follows, it 
is important to note that the interlayer bias also breaks 
inversion symmetry, and therefore, in combination with 
the intrinsic SOI, can lead to a spin splitting. 

The SOI in bilayer graphene is still a topic of on- 
going theoretical discussion [12, 13]. The Hamiltonian 
of the intrinsic SOI consistent with the crystal symme- 
try is found to be [12] H^o — Aitct 2 s z + \2T[i z s z + 
s x Tu x Sy) + \4& z (i-iyS x + T/jj X Sy), where /i^, <Ji 
and Si are Pauli matrices denoting layer, sublattice, and 
electron spin, respectively. The last SOI parameter which 
is estimated to be A4 = 0.48 meV dominates the other 
terms, with Ai = 14 \ieV, A2 = 8 ueV, and A3 = 5.5 |J.eV. 
Both the Ai and the A2 terms are diagonal in spin, pseu- 
dospin and layer leading to out-of-plane low-energy ef- 
fective spin-orbit fields which do not efficiently couple to 
momentum scattering as is needed for D'yakonov-Perel'- 
type spin relaxation. The remaining two terms give rise 
to in-plane spin-orbit fields which change their direc- 
tion depending on the angle of the electron's momentum. 
However, not only was A3 found to be much smaller than 
A4 in Ref. [12], but in comparison with A4-type spin-orbit 
interaction the magnitude of its corresponding spin-orbit 
field at low Fermi energies Ep is further supressed by 
Ep/ji. Below, we focus on the A4 term, for a discussion 
of the remaining terms of -ffso including the correspond- 
ing expressions for the spin-orbit fields, see Appendix. 

In the presence of SOI and for U 7^ 0, the four spin- 
degenerate bands described by Hblg split up into eight 
bands. Half of those bands are split off from the Dirac 
points by 71 and are not directly involved in spin trans- 
port when the Fermi energy is in the vicinity of the Dirac 
point. Among the remaining four low-energy bands, 
two correspond to electron and two to hole states, each 



2 



with their split spin degeneracy. To obtain the spin- 
orbit field for electrons (holes), we focus on positive 
(negative) Fermi energies, where spin currents are car- 
ried by the electrons (holes). In order to derive an ef- 
fective model for the low-energy bands, we perform a 
Schrieffcr- Wolff transformation on the total Hamiltonian 
H = HblgQ^S+Hso, restricting ourselves to the domi- 
nant A4 term for the rest of the discussion (see Appendix 
for a more general discussion). For this purpose, we di- 
vide up the total Hamiltonian into low- and high-energy 
parts (separated by 71), and the interactions V that cou- 
ple them, H = Hq + V, where H corresponds to Hblg 
without intralayer hopping (vf = 0), while V contains 
both intralayer hopping and SOI and can be expressed 
in the basis Ai it , A lti , B 2 ^, £2,4.1 B i,t> B i,U ^2,f, A 
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where & T ,±\ = (1± T )/2- 

We now perform the 
mation H = e s He~ s — 
anti-Hermitian matrix S 
the condition V + [S,H ] 
tions of order (\p\v F /-fi) 3 
neglected. The spin-independent part H° obtained 
from H by setting A^ = for i = 1,...,4 reproduces 
the known form of the low-energy bands [10]. E± = 



Schrieffer- Wolff transfor- 
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where k = hkyp/li and the (unnormalized) eigenstates 
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in the absence of SOI. The spin-dependent part H x = 
H — H° can be expressed in the eigenbasis Eq. (3) of H°, 
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with the electron (hole) effective spin-orbit field 
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The spin-orbit field and splitting are shown in Fig. 1 for 
two different values of the bias voltage, U = 0.1 eV and 
U = 0.01 cV. For A 4 < U, the SOI-induced electron-hole 
coupling A can be neglected, which is confirmed by our 
numerical analysis (see Fig. lb). 

As our next step, we derive the in- and out-of-plane 
spin relaxation times originating from the presence of 
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FIG. 1: (a) Spin-orbit field 17+ (k) of electrons in bilayer 
graphene. (b) Spin splitting AE where circles/diamonds 
refer to the energy difference between the two electron-like 
low-energy bands AE = |-E+,i — £+,2! obtained from a nu- 
merical diagonalization of the full Hamiltonian including A4- 
like SOI and lines to the splitting of the spin-orbit field 
AE = ft|ft+(k)|/2 given by Eq. (5). 



Ok = J~i + (k) via the D'yakonov-Perel' mechanism. For 
concreteness, we restrict ourselves to electrons. Spin 
transport is modeled using a kinetic spin Bloch equation 
(KSBE), i.e., a semiclassical rate equation for the spin 
distribution Sk carried by an ensemble of band electrons, 
an approach well known from semiconductor spintronics 
(see e.g. [14]). In the absence of external forces, 
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For the purpose of extracting the spin coherence times, it 
suffices to consider the simplified scenario of a homoge- 
neous spin distribution. Furthermore, we restrict our cal- 
culation to elastic, i.e., energy conserving, scattering and 
focus on the spin of the charge carriers at the Fermi sur- 
face, which essentially corresponds to a zero-temperature 
estimate. We consider Fermi energies Ep much smaller 
than energy separation of the split-off bands, but larger 
than the |k| = offset of the low energy bands, i.e. 
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71 3> £f > U/2. In this case there is a single connected 
Fermi surface near each of the two valleys K and K' and 
we can employ our effective low-energy theory with the 
spin-orbit field Eq. (5). 

At low energies, the energy bands experience a non- 
negligible anisotropy due to the trigonal warping intro- 
duced by the interlayer velocity V3, which substantially 
complicates solving the KSBE. However, the correspond- 
ing effect on spin-relaxation is in most cases relatively 
small which allows us to begin with a u 3 = estimate, 
subsequently include V3 to first order, and finally compare 
our analytical results to a numerical calculation taking 
trigonal warping fully into account. A description of the 
last two steps as well as a discussion of the different re- 
sults can be found in the Appendix. Here we only discuss 
the V3 = estimate. Deviations from this result are com- 
parably small and occur predominantly where the Fermi 
energy is low and very close to U/2. 

In the isotropic limit V3 — 0, the spin-orbit fields 
fi±(k) given for electrons in Eq. (5) and the band ener- 
gies E± simplify considerably. In particular, the magni- 
tude of the spin-orbit field becomes isotropic in this case, 
l^±(k)| = &±(k), and is therefore constant if we consider 
electrons at the Fermi level |p| — pf = hkp. Moreover, 
in this limit the spin-orbit field becomes independent of 
the valley. We can thus simply parameterize the spin dis- 
tribution in both valleys by the same angle <f>. In other 
words, it (formally) does not matter if the quasiparticle 
carrying the spin is located at K or K'. For elastic and 
symmetric scattering the scattering rates in Eq. (6) are of 
the form VF k ,k' = W k / jk = W(<f> - <j>')2-n:hv F 5(E k - E w ). 
In the isotropic limit the collision integral only needs to 
be taken over a circle of radius |k| = kp and the KSBE 
given by Eq. (6) reduces to 
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In order to solve Eq. (7) we first decompose the spin 
distribution function into an average So over the Fermi 
surface, which is independent of the angle </>, and the 
remaining deviation As k , describing the angular depen- 
dence, 
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s k = s + As k , 



so 



Sk 



2tt 



Sk, 



(8) 



where Note that the experimentally observed spin relax- 
ation refers to the decay of the total spin of the charge 
carriers at the Fermi surface, which is in turn given by 
the average spin polarization So. To obtain the time de- 
pendence of So we substitute Eq. (8) into the KSBE (7) 
and take the average over the angle cf>, 
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Note that both the spin-orbit field and the collision inte- 
gral average to zero. The corresponding equation for the 
anisotropic part is 
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The two coupled differential equations (9) and (10) can 
be solved approximately in the strong scattering limit 
I S~2 k I Tp <C 1, where r p denotes the momentum relaxation 
time. In this limit the combination of fast momentum 
scattering and slow spin precession implies that the de- 
viation Ask reaches a quasi-stationary state Asf when 
dAs^/dt w 0, which is then followed by a slow decay 
of the isotropic spin polarization Sq [14]. Since momen- 
tum relaxation is usually very fast on the time scale of 
the observation length, Ai D b s t p , the observed dynam- 
ics of the spin polarization s is effectively the averaged 
quantity sg bs (t) = J^J^t"/? s o{t')dt'. We can therefore 
neglect fast fluctuations occurring on the time scale r p 
as long as they are uncorrelated for times much longer 
than T p . It can be shown that the last two terms of 
Eq. (10) only give rise to fluctuations of the spin dis- 
tribution, which are uncorrelated on a time scale r p . 
Neglecting the last two terms of Eq. (10) the steady state 
condition becomes 

ttk x so = J ^W(<j> - </>') (Asf - Astf) . (11) 


This equation can be solved using the following ansatz, 

Asf = t* [Ok x s ] , (12) 

where we still need to determine the time r*. After 
substituting Eq. (12) into Eq. (11), the integral can be 
treated by expanding the scattering rates W((f> — </>') in 
polar harmonics. We find that Eq. (11) has the solution 
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and therefore r* can be identified with the momentum 
relaxation time, r* = r p . Having solved Eq. (11), we 
substitute the steady state solution Eq. (12) with Eq. (13) 
into the equation of motion of the total spin polarization 
Sq, Eq. (9), and find an exponential decay law, 
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FIG. 2: In-plane spin relaxation time Tg ii in bilayer graphene 
with A4-type SOI in the isotropic limit («3 = 0), (a) as a 
function of the Fermi energy Ef for different bias voltages U, 
and (b) as a function of the bias voltage U at a constant Fermi 
energy Ep = 0.06 eV. In both plots the momentum relaxation 
time is chosen to be r p = 10~~ 13 s. Dashed lines correspond to 
the large Fermi energy approximation given in Eq. (17). 
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is found by solving Ep — E+\ V3= q. The transverse spin 
relaxation time is simply ts,a_ = rsii/2. Combining 
Eqs. (15) and (16), we obtain the spin relaxation time 
as a function of the Fermi energy Ep and the bias volt- 
age U. As shown in Fig. 2, the spin relaxation time is 
very sensitive to both Ep and U . For a constant U and 
sufficiently large Ep the spin relaxation time increases as 
a function of Ep and can be approximated by [15] 
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The typical D'yakonov-Perel' relation 1/t$ oc t p has 
already been observed in two different experiments [9] 
and [8]. As pointed out in the previous discussion, the 
calculated relaxation rates are very sensitive to a number 



of parameters: the Fermi energy Ep, the bias voltage U, 
and the SOI strength A4. Unfortunately, these are not 
easily accessible experimentally. Thus, we will have to 
rely on some rough estimates in order to compare the 
experimental values of t$ with the obtained theoretical 
results. 

In [9], the spin relaxation time has been measured for 
for a range of mobilities from 300 cm 2 /Vs to 2000 cm 2 /Vs 
at room temperature and a range from 700cm 2 /Vs to 
3800 cm 2 /Vs at 5 K; both at a fixed carrier density n e — 
1.5T0 12 cm" 2 . We can roughly estimate the Fermi energy 
using the parabolic approximation E(k) m h 2 k 2 /2m* , 
where m* — ^\/2v 2 F is the effective mass [10]. Inte- 
grating the density of states D(E), which is constant 
within this approximation, one obtains a carrier den- 
sity of n e = m*Ep/TTh 2 . The estimated Fermi energy 
for the experimental carrier density n e = 1.5 • 10 12 cm~ 2 
is Ep = 67meV. Again, in the effective mass approxi- 
mation we can estimate the momentum relaxation time 
from t p = m* n/e [9]. Assuming that the bias offset U/2 
at k — is well below the Fermi energy (which is also 
necessary for the parabolic approximation) we can use 
the approximate spin relaxation rate given by equation 
(17), 
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The experimental results show a reasonable agreement 
with the model estimate for a bias voltage of 50meV. 
The corresponding model prediction Eq. (18) is T5 w 

n _ 1000 cm 2 /Vs 

0.5 ns ' — . 

At 5 K the bias voltage of 50 meV is significantly larger 
than the thermal energy (ksT s» 0.5 meV) and the Fermi 
energy of 67meV reasonably well above U/2. At room 
temperature (ksT sa 25 meV) the thermal energy is com- 
parable with both, thus making the zero temperature es- 
timate very approximate. Interestingly, the experimental 
data shows a stronger correlation t$ oh 1/fi at room tem- 
perature. 

In conclusion, we have calculated the spin relaxation 
time in bilayer graphene in dependence of the Fermi en- 
ergy and intcrlaycr bias potential. These two parameters 
can be tuned independently with top and back gates. 
Using experimentally determined parameters and mak- 
ing reasonable assumptions for the unknown values of 
U and Ep, we obtain good agreement with the existing 
experiments. We find a strong dependence of the spin re- 
laxation time on externally applied fields that may have 
applications in field-controlled spin valve devices. 
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Additional spin-orbit fields 

In the main text we have focused on A4-type SOI and derived the corresponding spin-orbit field. Analogous spin- 
orbit fields can however be derived for the omitted terms of Hgo, i.e. Ai-, A2-, and A3-type SOI. In lowest order of 
the spin-orbit coupling constants we can consider each of the above spin-orbit terms separately. For each term we can 
derive an analytic expression of the respective spin-orbit field using the same recipe as in the case of A4-type SOI. In 
order to compute (for i = l,2, 3) we start with H = Hblg ® Is + #so \\ s =o for j^i, i-e. we omit all terms in Hgo 
except for the one involving A;. Via a Schrieffer- Wolff transformation we separate high and low energy bands arriving 
at a effective low-energy Hamiltonian H, where we again neglect terms of order (p«f/7i) 3 or (Ai/71) 2 and higher. The 
resulting effective Hamiltonian H can be split into a kinetic and a spin-dependent part. In all three cases we recover 
the same spin- independent part part H° as previously for A4. The remaining spin-dependent part H x = H — H° 
is subsequently rotated into the Eigenbasis of H° as given by Eq. (3). For a sufficiently a large bias (U Ai), the 
electron-hole coupling A can be dropped. Form the remaining 2x2 blocks H*j h we obtain the respective spin-orbit 

field fJ±\ 

Below we report the resulting expressions for the approximate spin-orbit fields: 
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(21) 



Note that both fl^} and fl± are out-of-plane effective magnetic fields, which in the isotropic limit (^3 = 0) are 
independent of the electron momentum. Similar to fi^ 4 , fl^ is an in-plane effective field, which changes its direction 
depending on the angle of the electrons momentum (j>. In contrast to fi± 4 it is however not proportional to U jE\ , but 
instead to C//71. In other wor ds |0^ 3 |/|^ 4 | cx (A 3 /A 4 )(£;0/ 7l ) j w hich in the range of the low-energy theory makes 
it small even if A3 and A4 were comparable. 



Spin-relaxation - a first order estimate including trigonal warping 

In the analytic derivation of the in- and out-of-plane spin relaxation rates given in the main text we have neglected 
the anisotropy of the band structure. In the case of a finite trigonal warping ^ 0) the length of the Fermi 
wave vector is no longer constant on the Fermi surface. Moreover, the density of states at the Fermi level is no 
longer constant. Solving the general scattering integral, which previously used to be a simple integral over the angle, 
now becomes a more complicated task. In order to obtain a first estimate of the effect of trigonal warping on the 
spin relaxation time we instead choose a much simpler approach. Namely, we use the (momentum) relaxation time 
approximation of the KSBE. Here the scattering integral is replaced by a single parameter, the momentum relaxation 
time: 

— - fl k x s k = -t p (s k - (s k )) , (22) 

where (•) denotes the average over the Fermi surface and £7(k) = Ep , as we assume elastic scattering and only 
consider electrons at the Fermi level. Although we may not be able to calculate the average (•) analytically, we can 
use a semi-numerical approach. Therefore, we again decompose the spin distribution function into its average and its 
k-dependent deviation: 

s k = (s k ) + As k (23) 
Neglecting the fluctuation term Jl k x As k — (f2 k x As k ), the kinetic spin Bloch equation simplifies to 

%^ = (fikxAs k ) and ^ = ft k x (s k ) - r p As k . (24) 
ot at 
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The steady-state solution of the k-dependent part is readily given by As^' = t p (fit x (sk)). For the corresponding 
time evolution of the average spin we find 

= -t p (n* (s k ) st - (n k • ( Sk ) st ) n k ) . (25) 

The part of the right hand side that is proportional to (sk) st leads to the first order estimate of the spin relaxation 
time including trigonal warping, 

-j- = -r p (r! k -r! kil ) for i = x,y,z. (26) 

T S,i 

Since the out-of-plane component 17k, z simply vanishes and the two in-plane components are of the same average 
amplitude in both valleys, we can immediately recover that there is still only two different spin relaxation times (t^ m 
and rlJ. 

In order to numerically calculate Eq. (26), we need to explicitly calculate the average over the Fermi surface. In 
the case of 71 3> Ep > U/2 this can be achieved by a numerical inversion of the low-energy dispersion relation. 
Inverting E9_(k,<p) at a discrete number of angles and using a standard interpolating function, we obtain kp(tp), i.e. 
the amplitude of the Fermi vector as a function of its angle. The average over the Fermi surface can be expressed in 
terms of a single integral over the angle: 

</(M)> = \\ d<l>D(<t>)f[k F (ct>),cf>] , (27) 


where 

m= l V[M^ + [W ^ (28) 

is the respective density of states and Z = (D(4>)). The above density of states along the anisotropic Fermi surface 
can be derived form a coordinate transformation into local coordinates k» and k±, pointing along and perpendicular 
to the Fermi surface. 



A numerical model of spin relaxation 



To check the approximations we have employed when solving the KSBE, we also consider a simple numerical model 
that simulates the concept of the D'yakonov-Perel' mechanism. We therefore sample the spin evolution of an ensemble 
of electrons at the Fermi level. The diffusive (real space) motion of the ensemble is modeled by a random k-space 
walk of each electron. A homogeneous (or averaged) spin-orbit interaction is represented by a k-dependent spin-orbit 
field r2k, which in turn acts on the spin of each electron. Following the semiclassical approximation we assign each 
electron-like quasiparticle a wave vector k (relative to one of the Dirac-points) and a spin S. Their dynamics are 
governed by semiclassical equations of motion, i.e. in absence of external forces, unless the electron is being scattered, 
k is simply constant an S evolves according to dS/dt = f2 k x S. 

Momentum scattering on the other hand is modeled by a homogeneous scattering rate W(k, k'), representing the 
rate at which electrons in state k scatter into the state k'. Scattering is assumed to be elastic and spin conserving. 
For a simple model we consider scatterers to be represented by Gaussian model potentials of width R, i.e., V(r) = 
Vo exp(— r 2 /2R 2 ). Here we study small scatterers, where the spread of the potential is still larger than the lattice 
constant, but much smaller than the inverse of the wave vector amplitude: a <C R <C 1/k. In this limit the explicit 
|k|-dependence can be neglected and the scattering cross section simplifies to da/dO oc cos 2 9, where 9 is the scattering 
angle. The remaining dependence cos 2 9 is the signature of the Berry phase of the quasi particles. Note that R -C 1/k, 
where k <C K, also implies that intervalley scattering can be neglected. According to Fermi's golden rule the scattering 
rate form k to k' is proportional to the density of states at the outgoing momentum fck.' . If we again focus on Fermi 
energies sufficiently larger than the k = bias offset U/2, each electron wave vector k can be parameterized by its 
angle (f>, where |k| = kp(4>) (see previous section). This suggests the following expression for the scattering rate for 
small scatterers in bilayer graphene with finite trigonal warping, 

W(h<ft) = — \ cos 2 , (29) 



U [meV] E F [meV] r [ns] r\ [ns] r/« [ns] 

10 5l 3~li6 2^46 L99~~ 

25.1 4.26 4.25 4.20 

50.1 9.69 9.78 9.80 

50 25.1 0.547 0.623 0.612 

50.1 0.409 0.413 0.407 

100 50.1 0.171 0.188 0.190 



TABLE I: In-plane spin relaxation time in bilayer graphene with A4-type SOI for different bias voltages U and Fermi energies 
Ef- The table is a comparison of the zeroth and first order estimates, r§ n (15) and t§ y (26), with the best-fit values for the 
numerical data shown in Fig. 3. 



where Z = j d<fr'd<f) cos 2 (4> — <fi') D((f>') is the normalization, l/r sc the total scattering rate and D(4>') the angular 
dependent density of states as given by Eq. (28). All of the numerical results presented in Fig. 3 are calculated using 
this approximation. Note that in the isotropic limit D(4>') = const, implies that the momentum relaxation time t p 
(see Eq. 13) is equal to the mean scattering time r sc . 



Fermi energy and bias voltage dependence - a comparison with numerics 

Fig. 3 shows the numerical spin relaxation in comparison with the two estimates for a range of bias voltages from 
lOmeV to 100 meV and different Fermi energies. As previously noted, all Fermi energies are chosen to be larger than 
the k = offset given by the bias voltage (Ep > U/2). The calculated examples demonstrate an excellent agreement 
between the numerical data and the above first order estimate y (26) for all (b through f ) but the first example 
(a). In these cases even the zeroth order estimate T^y, where r p — t sc , is in comparably good agreement with the 
numerical data. Noteworthy deviations only occur in the cases (d and f), where the Fermi energy is very close to the 
voltage offset (Ep — U/2 + 0.1 meV). As shown in the corresponding insets these are exactly the cases where trigonal 
warping is most pronounced, (a) is the only example where both t§ y and Tg y deviate significantly from the numerical 
results. However, taking into account the extreme trigonal warping, both still provide a good order of magnitude 
estimate. Overall, the numerical data supports the sensitive dependence of the spin relaxation time on bias voltage 
and Fermi energy shown in Fig. 2. Notice that there is roughly two orders of magnitude difference between the spin 
relaxation times for U = lOmeV and U — 100 meV at Ep = 50.1meV (see table I). 

Though not explicitly shown here, the numerical calculations indicate the same anisotropy factor of two between 
the relaxation times of the in- and out-of-plane spin polarization, which we derived analytically in the isotropic limit 

V 3 = 0. , 
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a) U = lOmeV E F = 5.1 meV d) U = 50 meV E F = 25.1 meV 




1 2 3 4 5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 

t [ns] t [ns] 



b) U = lOmeV E F = 25.1 meV e) U = 50 meV E F = 50.1 meV 




1 2 3 4 5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 

t [ns] t [ns] 



c) £/ = lOmeV E F = 50.1 meV f) U= 100 meV E F = 50.1 meV 




1 2 3 4 5 0.0 0.1 0.2 0.3 0.4 0.5 

f [ns] f [ns] 



FIG. 3: (color online) In-plane spin relaxation in bilayer graphene with A4-type SOI. The plots show a comparison of the 
numerical data (diamonds) and the zeroth (orange dashed) and first (continuous purple lines) order estimates, m (15) and 
Tg || (26), for different values of the interlayer bias U and the Fermi energy E F . All curves are calculated in the limit of 
a «C R <g l/k F , using a mean scattering time r sc = r p = 0.1 ps. The Insets show the trigonal warping of the Fermi surfaces. 
The corresponding spin relaxation times are listed in table I. 



